home *** CD-ROM | disk | FTP | other *** search
/ Aminet 30 / Aminet 30 (1999)(Schatztruhe)[!][Apr 1999].iso / Aminet / gfx / misc / gnuplot-3.7src.lha / gnuplot-3.7src / gnuplot-3.7.lha / gnuplot-3.7 / fit.c < prev    next >
C/C++ Source or Header  |  1999-01-06  |  45KB  |  1,609 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: fit.c,v 1.58 1998/04/14 00:15:19 drd Exp $";
  3. #endif
  4.  
  5. /*
  6.  *    Nonlinear least squares fit according to the
  7.  *      Marquardt-Levenberg-algorithm
  8.  *
  9.  *      added as Patch to Gnuplot (v3.2 and higher)
  10.  *      by Carsten Grammes
  11.  *      Experimental Physics, University of Saarbruecken, Germany
  12.  *
  13.  *      Internet address: cagr@rz.uni-sb.de
  14.  *
  15.  *      Copyright of this module:  1993, 1998  Carsten Grammes
  16.  *
  17.  *      Permission to use, copy, and distribute this software and its
  18.  *      documentation for any purpose with or without fee is hereby granted,
  19.  *      provided that the above copyright notice appear in all copies and
  20.  *      that both that copyright notice and this permission notice appear
  21.  *      in supporting documentation.
  22.  *
  23.  *      This software is provided "as is" without express or implied warranty.
  24.  *
  25.  *      930726:     Recoding of the Unix-like raw console I/O routines by:
  26.  *                  Michele Marziani (marziani@ferrara.infn.it)
  27.  * drd: start unitialised variables at 1 rather than NEARLY_ZERO
  28.  *  (fit is more likely to converge if started from 1 than 1e-30 ?)
  29.  *
  30.  * HBB (Broeker@physik.rwth-aachen.de) : fit didn't calculate the errors
  31.  * in the 'physically correct' (:-) way, if a third data column containing
  32.  * the errors (or 'uncertainties') of the input data was given. I think
  33.  * I've fixed that, but I'm not sure I really understood the M-L-algo well
  34.  * enough to get it right. I deduced my change from the final steps of the
  35.  * equivalent algorithm for the linear case, which is much easier to
  36.  * understand. (I also made some minor, mostly cosmetic changes)
  37.  *
  38.  * HBB (again): added error checking for negative covar[i][i] values and
  39.  * for too many parameters being specified.
  40.  *
  41.  * drd: allow 3d fitting. Data value is now called fit_z internally,
  42.  * ie a 2d fit is z vs x, and a 3d fit is z vs x and y.
  43.  *
  44.  * Lars Hecking : review update command, for VMS in particular, where
  45.  * it is not necessary to rename the old file.
  46.  *
  47.  * HBB, 971023: lifted fixed limit on number of datapoints, and number
  48.  * of parameters.
  49.  */
  50.  
  51.  
  52. #define FIT_MAIN
  53.  
  54. #include <signal.h>
  55.  
  56. #include "plot.h"
  57. #include "matrix.h"
  58. #include "fit.h"
  59. #include "setshow.h"        /* for load_range */
  60. #include "alloc.h"
  61.  
  62. #if defined(MSDOS) || defined(DOS386)    /* non-blocking IO stuff */
  63. # include <io.h>
  64. # ifndef _Windows        /* WIN16 does define MSDOS .... */
  65. #  include <conio.h>
  66. # endif
  67. # include <dos.h>
  68. #else /* !(MSDOS || DOS386) */
  69. # ifndef VMS
  70. #  include <fcntl.h>
  71. # endif                /* !VMS */
  72. #endif /* !(MSDOS || DOS386) */
  73.  
  74. #if defined(ATARI) || defined(MTOS)
  75. # define getchx() Crawcin()
  76. static int kbhit(void);
  77. #endif
  78.  
  79. #define STANDARD    stderr    /* compatible with gnuplot philosophy */
  80.  
  81. #define BACKUP_SUFFIX ".old"
  82.  
  83.  
  84. /* access external global variables  (ought to make a globals.h someday) */
  85.  
  86. extern struct udft_entry *dummy_func;
  87. extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  88. extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  89. extern int c_token;
  90. extern int df_datum, df_line_number;
  91.  
  92.  
  93. enum marq_res {
  94.     OK, ERROR, BETTER, WORSE
  95. };
  96. typedef enum marq_res marq_res_t;
  97.  
  98. #ifdef INFINITY
  99. # undef INFINITY
  100. #endif
  101.  
  102. #define INFINITY    1e30
  103. #define NEARLY_ZERO 1e-30
  104.  
  105. /* create new variables with this value (was NEARLY_ZERO) */
  106. #define INITIAL_VALUE 1.0
  107.  
  108. /* Relative change for derivatives */
  109. #define DELTA        0.001
  110.  
  111. #define MAX_DATA    2048
  112. #define MAX_PARAMS  32
  113. #define MAX_LAMBDA  1e20
  114. #define MIN_LAMBDA  1e-20
  115. #define LAMBDA_UP_FACTOR 10
  116. #define LAMBDA_DOWN_FACTOR 10
  117.  
  118. #if defined(MSDOS) || defined(OS2) || defined(DOS386)
  119. # define PLUSMINUS   "\xF1"    /* plusminus sign */
  120. #else
  121. # define PLUSMINUS   "+/-"
  122. #endif
  123.  
  124. #define LASTFITCMDLENGTH 511
  125.  
  126. /* HBB 971023: new, allow for dynamic adjustment of these: */
  127. static int max_data;
  128. static int max_params;
  129.  
  130. static double epsilon = 1e-5;    /* convergence limit */
  131. static int maxiter = 0;        /* HBB 970304: maxiter patch */
  132.  
  133. static char fit_script[128];
  134. static char logfile[128] = "fit.log";
  135. static char *FIXED = "# FIXED";
  136. static char *GNUFITLOG = "FIT_LOG";
  137. static char *FITLIMIT = "FIT_LIMIT";
  138. static char *FITSTARTLAMBDA = "FIT_START_LAMBDA";
  139. static char *FITLAMBDAFACTOR = "FIT_LAMBDA_FACTOR";
  140. static char *FITMAXITER = "FIT_MAXITER";    /* HBB 970304: maxiter patch */
  141. static char *FITSCRIPT = "FIT_SCRIPT";
  142. static char *DEFAULT_CMD = "replot";    /* if no fitscript spec. */
  143. static char last_fit_command[LASTFITCMDLENGTH+1] = "";
  144.  
  145.  
  146. static FILE *log_f = NULL;
  147.  
  148. static int num_data, num_params;
  149. static int columns;
  150. static double *fit_x = 0, *fit_y = 0, *fit_z = 0, *err_data = 0, *a = 0;
  151. static TBOOLEAN ctrlc_flag = FALSE;
  152. static TBOOLEAN user_stop = FALSE;
  153.  
  154. static struct udft_entry func;
  155.  
  156. typedef char fixstr[MAX_ID_LEN+1];
  157.  
  158. static fixstr *par_name;
  159.  
  160. static double startup_lambda = 0, lambda_down_factor = LAMBDA_DOWN_FACTOR,
  161.  lambda_up_factor = LAMBDA_UP_FACTOR;
  162.  
  163. /*****************************************************************
  164.              internal Prototypes
  165. *****************************************************************/
  166.  
  167. static RETSIGTYPE ctrlc_handle __PROTO((int an_int));
  168. static void ctrlc_setup __PROTO((void));
  169. static void printmatrix __PROTO((double **C, int m, int n));
  170. static void print_matrix_and_vectors __PROTO((double **C, double *d, double *r, int m, int n));
  171. static marq_res_t marquardt __PROTO((double a[], double **alpha, double *chisq,
  172.                      double *lambda));
  173. static TBOOLEAN analyze __PROTO((double a[], double **alpha, double beta[],
  174.                  double *chisq));
  175. static void calculate __PROTO((double *zfunc, double **dzda, double a[]));
  176. static void call_gnuplot __PROTO((double *par, double *data));
  177. static TBOOLEAN fit_interrupt __PROTO((void));
  178. static TBOOLEAN regress __PROTO((double a[]));
  179. static void show_fit __PROTO((int i, double chisq, double last_chisq, double *a,
  180.                   double lambda, FILE * device));
  181. static TBOOLEAN is_empty __PROTO((char *s));
  182. static TBOOLEAN is_variable __PROTO((char *s));
  183. static double getdvar __PROTO((char *varname));
  184. static double createdvar __PROTO((char *varname, double value));
  185. static void splitpath __PROTO((char *s, char *p, char *f));
  186. static void backup_file __PROTO((char *, const char *));
  187.  
  188.  
  189. /*****************************************************************
  190.     Small function to write the last fit command into a file
  191.     Arg: Pointer to the file; if NULL, nothing is written,
  192.          but only the size of the string is returned.
  193. *****************************************************************/
  194.  
  195. size_t wri_to_fil_last_fit_cmd(fp)
  196. FILE *fp;
  197. {
  198.     if (fp == NULL)
  199.     return strlen(last_fit_command);
  200.     else
  201.     return (size_t) fputs(last_fit_command, fp);
  202. }
  203.  
  204.  
  205. /*****************************************************************
  206.     This is called when a SIGINT occurs during fit
  207. *****************************************************************/
  208.  
  209. static RETSIGTYPE ctrlc_handle(an_int)
  210. int an_int;
  211. {
  212. #ifdef OS2
  213.     (void) signal(an_int, SIG_ACK);
  214. #else
  215.     /* reinstall signal handler (necessary on SysV) */
  216.     (void) signal(SIGINT, (sigfunc) ctrlc_handle);
  217. #endif
  218.     ctrlc_flag = TRUE;
  219. }
  220.  
  221.  
  222. /*****************************************************************
  223.     setup the ctrl_c signal handler
  224. *****************************************************************/
  225. static void ctrlc_setup()
  226. {
  227. /*
  228.  *  MSDOS defines signal(SIGINT) but doesn't handle it through
  229.  *  real interrupts. So there remain cases in which a ctrl-c may
  230.  *  be uncaught by signal. We must use kbhit() instead that really
  231.  *  serves the keyboard interrupt (or write an own interrupt func
  232.  *  which also generates #ifdefs)
  233.  *
  234.  *  I hope that other OSes do it better, if not... add #ifdefs :-(
  235.  */
  236. #if (defined(__EMX__) || !defined(MSDOS) && !defined(DOS386)) && !defined(ATARI) && !defined(MTOS)
  237.     (void) signal(SIGINT, (sigfunc) ctrlc_handle);
  238. #endif
  239. }
  240.  
  241.  
  242. /*****************************************************************
  243.     getch that handles also function keys etc.
  244. *****************************************************************/
  245. #if defined(MSDOS) || defined(DOS386)
  246.  
  247. /* HBB 980317: added a prototype... */
  248. int getchx __PROTO((void));
  249.  
  250. int getchx()
  251. {
  252.     register int c = getch();
  253.     if (!c || c == 0xE0) {
  254.     c <<= 8;
  255.     c |= getch();
  256.     }
  257.     return c;
  258. }
  259. #endif
  260.  
  261. /*****************************************************************
  262.     in case of fatal errors
  263. *****************************************************************/
  264. void error_ex()
  265. {
  266.     char *sp;
  267.  
  268.     strncpy(fitbuf, "         ", 9);        /* start after GNUPLOT> */
  269.     sp = strchr(fitbuf, NUL);
  270.     while (*--sp == '\n')
  271.     ;
  272.     strcpy(sp + 1, "\n\n");    /* terminate with exactly 2 newlines */
  273.     fputs(fitbuf, STANDARD);
  274.     if (log_f) {
  275.     fprintf(log_f, "BREAK: %s", fitbuf);
  276.     (void) fclose(log_f);
  277.     log_f = NULL;
  278.     }
  279.     if (func.at) {
  280.     free(func.at);        /* release perm. action table */
  281.     func.at = (struct at_type *) NULL;
  282.     }
  283.     /* restore original SIGINT function */
  284. #if !defined(ATARI) && !defined(MTOS)
  285.     interrupt_setup();
  286. #endif
  287.  
  288.     bail_to_command_line();
  289. }
  290.  
  291.  
  292. /*****************************************************************
  293.     New utility routine: print a matrix (for debugging the alg.)
  294. *****************************************************************/
  295. static void printmatrix(C, m, n)
  296. double **C;
  297. int m, n;
  298. {
  299.     int i, j;
  300.  
  301.     for (i = 0; i < m; i++) {
  302.     for (j = 0; j < n - 1; j++)
  303.         Dblf2("%.8g |", C[i][j]);
  304.     Dblf2("%.8g\n", C[i][j]);
  305.     }
  306.     Dblf("\n");
  307. }
  308.  
  309. /**************************************************************************
  310.     Yet another debugging aid: print matrix, with diff. and residue vector
  311. **************************************************************************/
  312. static void print_matrix_and_vectors(C, d, r, m, n)
  313. double **C;
  314. double *d, *r;
  315. int m, n;
  316. {
  317.     int i, j;
  318.  
  319.     for (i = 0; i < m; i++) {
  320.     for (j = 0; j < n; j++)
  321.         Dblf2("%8g ", C[i][j]);
  322.     Dblf3("| %8g | %8g\n", d[i], r[i]);
  323.     }
  324.     Dblf("\n");
  325. }
  326.  
  327.  
  328. /*****************************************************************
  329.     Marquardt's nonlinear least squares fit
  330. *****************************************************************/
  331. static marq_res_t marquardt(a, C, chisq, lambda)
  332. double a[];
  333. double **C;
  334. double *chisq;
  335. double *lambda;
  336. {
  337.     int i, j;
  338.     static double *da = 0,    /* delta-step of the parameter */
  339.     *temp_a = 0,        /* temptative new params set   */
  340.     *d = 0, *tmp_d = 0, **tmp_C = 0, *residues = 0;
  341.     double tmp_chisq;
  342.  
  343.     /* Initialization when lambda == -1 */
  344.  
  345.     if (*lambda == -1) {    /* Get first chi-square check */
  346.     TBOOLEAN analyze_ret;
  347.  
  348.     temp_a = vec(num_params);
  349.     d = vec(num_data + num_params);
  350.     tmp_d = vec(num_data + num_params);
  351.     da = vec(num_params);
  352.     residues = vec(num_data + num_params);
  353.     tmp_C = matr(num_data + num_params, num_params);
  354.  
  355.     analyze_ret = analyze(a, C, d, chisq);
  356.  
  357.     /* Calculate a useful startup value for lambda, as given by Schwarz */
  358.     /* FIXME: this is doesn't turn out to be much better, really... */
  359.     if (startup_lambda != 0)
  360.         *lambda = startup_lambda;
  361.     else {
  362.         *lambda = 0;
  363.         for (i = 0; i < num_data; i++)
  364.         for (j = 0; j < num_params; j++)
  365.             *lambda += C[i][j] * C[i][j];
  366.         *lambda = sqrt(*lambda / num_data / num_params);
  367.     }
  368.  
  369.     /* Fill in the lower square part of C (the diagonal is filled in on
  370.        each iteration, see below) */
  371.     for (i = 0; i < num_params; i++)
  372.         for (j = 0; j < i; j++)
  373.         C[num_data + i][j] = 0, C[num_data + j][i] = 0;
  374.     /*printmatrix(C, num_data+num_params, num_params); */
  375.     return analyze_ret ? OK : ERROR;
  376.     }
  377.     /* once converged, free dynamic allocated vars */
  378.  
  379.     if (*lambda == -2) {
  380.     free(d);
  381.     free(tmp_d);
  382.     free(da);
  383.     free(temp_a);
  384.     free(residues);
  385.     free_matr(tmp_C);
  386.     return OK;
  387.     }
  388.     /* Givens calculates in-place, so make working copies of C and d */
  389.  
  390.     for (j = 0; j < num_data + num_params; j++)
  391.     memcpy(tmp_C[j], C[j], num_params * sizeof(double));
  392.     memcpy(tmp_d, d, num_data * sizeof(double));
  393.  
  394.     /* fill in additional parts of tmp_C, tmp_d */
  395.  
  396.     for (i = 0; i < num_params; i++) {
  397.     /* fill in low diag. of tmp_C ... */
  398.     tmp_C[num_data + i][i] = *lambda;
  399.     /* ... and low part of tmp_d */
  400.     tmp_d[num_data + i] = 0;
  401.     }
  402.     /* printmatrix(tmp_C, num_data+num_params, num_params); */
  403.  
  404.     /* FIXME: residues[] isn't used at all. Why? Should it be used? */
  405.  
  406.     Givens(tmp_C, tmp_d, da, residues, num_params + num_data, num_params, 1);
  407.     /*print_matrix_and_vectors (tmp_C, tmp_d, residues,
  408.        num_params+num_data, num_params); */
  409.  
  410.     /* check if trial did ameliorate sum of squares */
  411.  
  412.     for (j = 0; j < num_params; j++)
  413.     temp_a[j] = a[j] + da[j];
  414.  
  415.     if (!analyze(temp_a, tmp_C, tmp_d, &tmp_chisq)) {
  416.     /* FIXME: will never be reached: always returns TRUE */
  417.     return ERROR;
  418.     }
  419.  
  420.     if (tmp_chisq < *chisq) {    /* Success, accept new solution */
  421.     if (*lambda > MIN_LAMBDA) {
  422.         (void) putc('/', stderr);
  423.         *lambda /= lambda_down_factor;
  424.     }
  425.     *chisq = tmp_chisq;
  426.     for (j = 0; j < num_data; j++) {
  427.         memcpy(C[j], tmp_C[j], num_params * sizeof(double));
  428.         d[j] = tmp_d[j];
  429.     }
  430.     for (j = 0; j < num_params; j++)
  431.         a[j] = temp_a[j];
  432.     return BETTER;
  433.     } else {            /* failure, increase lambda and return */
  434.     (void) putc('*', stderr);
  435.     *lambda *= lambda_up_factor;
  436.     return WORSE;
  437.     }
  438. }
  439.  
  440.  
  441. /* FIXME: in the new code, this function doesn't really do enough to be
  442.  * useful. Maybe it ought to be deleted, i.e. integrated with
  443.  * calculate() ?
  444.  */
  445. /*****************************************************************
  446.     compute chi-square and numeric derivations
  447. *****************************************************************/
  448. static TBOOLEAN analyze(a, C, d, chisq)
  449. double a[];
  450. double **C;
  451. double d[];
  452. double *chisq;
  453. {
  454. /*
  455.  *  used by marquardt to evaluate the linearized fitting matrix C
  456.  *  and vector d, fills in only the top part of C and d
  457.  *  I don't use a temporary array zfunc[] any more. Just use
  458.  *  d[] instead.
  459.  */
  460.     int i, j;
  461.  
  462.     *chisq = 0;
  463.     calculate(d, C, a);
  464.  
  465.     for (i = 0; i < num_data; i++) {
  466.     /* note: order reversed, as used by Schwarz */
  467.     d[i] = (d[i] - fit_z[i]) / err_data[i];
  468.     *chisq += d[i] * d[i];
  469.     for (j = 0; j < num_params; j++)
  470.         C[i][j] /= err_data[i];
  471.     }
  472.     /* FIXME: why return a value that is always TRUE ? */
  473.     return TRUE;
  474. }
  475.  
  476.  
  477. /* To use the more exact, but slower two-side formula, activate the
  478.    following line: */
  479. /*#define TWO_SIDE_DIFFERENTIATION */
  480. /*****************************************************************
  481.     compute function values and partial derivatives of chi-square
  482. *****************************************************************/
  483. static void calculate(zfunc, dzda, a)
  484. double *zfunc;
  485. double **dzda;
  486. double a[];
  487. {
  488.     int k, p;
  489.     double tmp_a;
  490.     double *tmp_high, *tmp_pars;
  491. #ifdef TWO_SIDE_DIFFERENTIATION
  492.     double *tmp_low;
  493. #endif
  494.  
  495.     tmp_high = vec(num_data);    /* numeric derivations */
  496. #ifdef TWO_SIDE_DIFFERENTIATION
  497.     tmp_low = vec(num_data);
  498. #endif
  499.     tmp_pars = vec(num_params);
  500.  
  501.     /* first function values */
  502.  
  503.     call_gnuplot(a, zfunc);
  504.  
  505.     /* then derivatives */
  506.  
  507.     for (p = 0; p < num_params; p++)
  508.     tmp_pars[p] = a[p];
  509.     for (p = 0; p < num_params; p++) {
  510.     tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
  511.     tmp_pars[p] = tmp_a * (1 + DELTA);
  512.     call_gnuplot(tmp_pars, tmp_high);
  513. #ifdef TWO_SIDE_DIFFERENTIATION
  514.     tmp_pars[p] = tmp_a * (1 - DELTA);
  515.     call_gnuplot(tmp_pars, tmp_low);
  516. #endif
  517.     for (k = 0; k < num_data; k++)
  518. #ifdef TWO_SIDE_DIFFERENTIATION
  519.         dzda[k][p] = (tmp_high[k] - tmp_low[k]) / (2 * tmp_a * DELTA);
  520. #else
  521.         dzda[k][p] = (tmp_high[k] - zfunc[k]) / (tmp_a * DELTA);
  522. #endif
  523.     tmp_pars[p] = a[p];
  524.     }
  525.  
  526. #ifdef TWO_SIDE_DIFFERENTIATION
  527.     free(tmp_low);
  528. #endif
  529.     free(tmp_high);
  530.     free(tmp_pars);
  531. }
  532.  
  533.  
  534. /*****************************************************************
  535.     call internal gnuplot functions
  536. *****************************************************************/
  537. static void call_gnuplot(par, data)
  538. double *par;
  539. double *data;
  540. {
  541.     int i;
  542.     struct value v;
  543.  
  544.     /* set parameters first */
  545.     for (i = 0; i < num_params; i++) {
  546.     (void) Gcomplex(&v, par[i], 0.0);
  547.     setvar(par_name[i], v);
  548.     }
  549.  
  550.     for (i = 0; i < num_data; i++) {
  551.     /* calculate fit-function value */
  552.     (void) Gcomplex(&func.dummy_values[0], fit_x[i], 0.0);
  553.     (void) Gcomplex(&func.dummy_values[1], fit_y[i], 0.0);
  554.     evaluate_at(func.at, &v);
  555.     if (undefined)
  556.         Eex("Undefined value during function evaluation");
  557.     data[i] = real(&v);
  558.     }
  559. }
  560.  
  561.  
  562. /*****************************************************************
  563.     handle user interrupts during fit
  564. *****************************************************************/
  565. static TBOOLEAN fit_interrupt()
  566. {
  567.     while (TRUE) {
  568.     fputs("\n\n(S)top fit, (C)ontinue, (E)xecute FIT_SCRIPT:  ", STANDARD);
  569.     switch (getc(stdin)) {
  570.  
  571.     case EOF:
  572.     case 's':
  573.     case 'S':
  574.         fputs("Stop.", STANDARD);
  575.         user_stop = TRUE;
  576.         return FALSE;
  577.  
  578.     case 'c':
  579.     case 'C':
  580.         fputs("Continue.", STANDARD);
  581.         return TRUE;
  582.  
  583.     case 'e':
  584.     case 'E':{
  585.         int i;
  586.         struct value v;
  587.         char *tmp;
  588.  
  589.         tmp = (fit_script != 0 && *fit_script) ? fit_script : DEFAULT_CMD;
  590.         fprintf(STANDARD, "executing: %s", tmp);
  591.         /* set parameters visible to gnuplot */
  592.         for (i = 0; i < num_params; i++) {
  593.             (void) Gcomplex(&v, a[i], 0.0);
  594.             setvar(par_name[i], v);
  595.         }
  596.         sprintf(input_line, tmp);
  597.         (void) do_line();
  598.         }
  599.     }
  600.     }
  601.     return TRUE;
  602. }
  603.  
  604.  
  605. /*****************************************************************
  606.     frame routine for the marquardt-fit
  607. *****************************************************************/
  608. static TBOOLEAN regress(a)
  609. double a[];
  610. {
  611.     double **covar, *dpar, **C, chisq, last_chisq, lambda;
  612.     int iter, i, j;
  613.     marq_res_t res;
  614.  
  615.     chisq = last_chisq = INFINITY;
  616.     C = matr(num_data + num_params, num_params);
  617.     lambda = -1;        /* use sign as flag */
  618.     iter = 0;            /* iteration counter  */
  619.  
  620.     /* ctrlc now serves as Hotkey */
  621.     ctrlc_setup();
  622.  
  623.     /* Initialize internal variables and 1st chi-square check */
  624.     if ((res = marquardt(a, C, &chisq, &lambda)) == ERROR)
  625.     Eex("FIT: error occured during fit");
  626.     res = BETTER;
  627.  
  628.     show_fit(iter, chisq, chisq, a, lambda, STANDARD);
  629.     show_fit(iter, chisq, chisq, a, lambda, log_f);
  630.  
  631.     /* MAIN FIT LOOP: do the regression iteration */
  632.  
  633.     /* HBB 981118: initialize new variable 'user_break' */
  634.     user_stop = FALSE;
  635.  
  636.     do {
  637. /*
  638.  *  MSDOS defines signal(SIGINT) but doesn't handle it through
  639.  *  real interrupts. So there remain cases in which a ctrl-c may
  640.  *  be uncaught by signal. We must use kbhit() instead that really
  641.  *  serves the keyboard interrupt (or write an own interrupt func
  642.  *  which also generates #ifdefs)
  643.  *
  644.  *  I hope that other OSes do it better, if not... add #ifdefs :-(
  645.  *  EMX does not have kbhit.
  646.  * 
  647.  *  HBB: I think this can be enabled for DJGPP V2. SIGINT is actually
  648.  *  handled there, AFAIK.
  649.  */
  650. #if ((defined(MSDOS) || defined(DOS386)) && !defined(__EMX__)) || defined(ATARI) || defined(MTOS)
  651.     if (kbhit()) {
  652.         do {
  653.         getchx();
  654.         } while (kbhit());
  655.         ctrlc_flag = TRUE;
  656.     }
  657. #endif
  658.  
  659.     if (ctrlc_flag) {
  660.         show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
  661.         ctrlc_flag = FALSE;
  662.         if (!fit_interrupt())    /* handle keys */
  663.         break;
  664.     }
  665.     if (res == BETTER) {
  666.         iter++;
  667.         last_chisq = chisq;
  668.     }
  669.     if ((res = marquardt(a, C, &chisq, &lambda)) == BETTER)
  670.         show_fit(iter, chisq, last_chisq, a, lambda, STANDARD);
  671.     } while ((res != ERROR)
  672.          && (lambda < MAX_LAMBDA)
  673.          && ((maxiter == 0) || (iter <= maxiter))
  674.          && (res == WORSE
  675.          || ((chisq > NEARLY_ZERO)
  676.              ? ((last_chisq - chisq) / chisq)
  677.              : (last_chisq - chisq)) > epsilon
  678.          )
  679.     );
  680.  
  681.     /* fit done */
  682.  
  683. #if !defined(ATARI) && !defined(MTOS)
  684.     /* restore original SIGINT function */
  685.     interrupt_setup();
  686. #endif
  687.  
  688.     /* HBB 970304: the maxiter patch: */
  689.     if ((maxiter > 0) && (iter > maxiter)) {
  690.     Dblf2("\nMaximum iteration count (%d) reached. Fit stopped.\n", maxiter);
  691.     } else if (user_stop) {
  692.     Dblf2("\nThe fit was stopped by the user after %d iterations.\n", iter);
  693.     } else {
  694.     Dblf2("\nAfter %d iterations the fit converged.\n", iter);
  695.     }
  696.  
  697.     Dblf2("final sum of squares of residuals : %g\n", chisq);
  698.     if (chisq > NEARLY_ZERO) {
  699.     Dblf2("rel. change during last iteration : %g\n\n", (chisq - last_chisq) / chisq);
  700.     } else {
  701.     Dblf2("abs. change during last iteration : %g\n\n", (chisq - last_chisq));
  702.     }
  703.  
  704.     if (res == ERROR)
  705.     Eex("FIT: error occured during fit");
  706.  
  707.     /* compute errors in the parameters */
  708.  
  709.     if (num_data == num_params) {
  710.     int i;
  711.  
  712.     Dblf("\nExactly as many data points as there are parameters.\n");
  713.     Dblf("In this degenerate case, all errors are zero by definition.\n\n");
  714.     Dblf("Final set of parameters \n");
  715.     Dblf("======================= \n\n");
  716.     for (i = 0; i < num_params; i++)
  717.         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]);
  718.     } else if (chisq < NEARLY_ZERO) {
  719.     int i;
  720.  
  721.     Dblf("\nHmmmm.... Sum of squared residuals is zero. Can't compute errors.\n\n");
  722.     Dblf("Final set of parameters \n");
  723.     Dblf("======================= \n\n");
  724.     for (i = 0; i < num_params; i++)
  725.         Dblf3("%-15.15s = %-15g\n", par_name[i], a[i]);
  726.     } else {
  727.         Dblf2("degrees of freedom (ndf) : %d\n",  num_data - num_params);
  728.         Dblf2("rms of residuals      (stdfit) = sqrt(WSSR/ndf)      : %g\n", sqrt(chisq / (num_data - num_params)));
  729.          Dblf2("variance of residuals (reduced chisquare) = WSSR/ndf : %g\n\n", chisq / (num_data - num_params));
  730.  
  731.     /* get covariance-, Korrelations- and Kurvature-Matrix */
  732.     /* and errors in the parameters                     */
  733.  
  734.     /* compute covar[][] directly from C */
  735.     Givens(C, 0, 0, 0, num_data, num_params, 0);
  736.     /*printmatrix(C, num_params, num_params); */
  737.  
  738.     /* Use lower square of C for covar */
  739.     covar = C + num_data;
  740.     Invert_RtR(C, covar, num_params);
  741.     /*printmatrix(covar, num_params, num_params); */
  742.  
  743.     /* calculate unscaled parameter errors in dpar[]: */
  744.     dpar = vec(num_params);
  745.     for (i = 0; i < num_params; i++) {
  746.         /* FIXME: can this still happen ? */
  747.         if (covar[i][i] <= 0.0)    /* HBB: prevent floating point exception later on */
  748.         Eex("Calculation error: non-positive diagonal element in covar. matrix");
  749.         dpar[i] = sqrt(covar[i][i]);
  750.     }
  751.  
  752.     /* transform covariances into correlations */
  753.     for (i = 0; i < num_params; i++) {
  754.         /* only lower triangle needs to be handled */
  755.         for (j = 0; j <= i; j++)
  756.         covar[i][j] /= dpar[i] * dpar[j];
  757.     }
  758.  
  759.     /* scale parameter errors based on chisq */
  760.     chisq = sqrt(chisq / (num_data - num_params));
  761.     for (i = 0; i < num_params; i++)
  762.         dpar[i] *= chisq;
  763.  
  764.     Dblf("Final set of parameters            Asymptotic Standard Error\n");
  765.     Dblf("=======================            ==========================\n\n");
  766.  
  767.     for (i = 0; i < num_params; i++) {
  768.         double temp =
  769.         (fabs(a[i]) < NEARLY_ZERO) ? 0.0 : fabs(100.0 * dpar[i] / a[i]);
  770.         Dblf6("%-15.15s = %-15g  %-3.3s %-12.4g (%.4g%%)\n",
  771.           par_name[i], a[i], PLUSMINUS, dpar[i], temp);
  772.     }
  773.  
  774.     Dblf("\n\ncorrelation matrix of the fit parameters:\n\n");
  775.     Dblf("               ");
  776.  
  777.     for (j = 0; j < num_params; j++)
  778.         Dblf2("%-6.6s ", par_name[j]);
  779.  
  780.     Dblf("\n");
  781.     for (i = 0; i < num_params; i++) {
  782.         Dblf2("%-15.15s", par_name[i]);
  783.         for (j = 0; j <= i; j++) {
  784.         /* Only print lower triangle of symmetric matrix */
  785.         Dblf2("%6.3f ", covar[i][j]);
  786.         }
  787.         Dblf("\n");
  788.     }
  789.  
  790.     free(dpar);
  791.     }
  792.  
  793.     /* call destructor for allocated vars */
  794.     lambda = -2;        /* flag value, meaning 'destruct!' */
  795.     (void) marquardt(a, C, &chisq, &lambda);
  796.  
  797.     free_matr(C);
  798.     return TRUE;
  799. }
  800.  
  801.  
  802. /*****************************************************************
  803.     display actual state of the fit
  804. *****************************************************************/
  805. static void show_fit(i, chisq, last_chisq, a, lambda, device)
  806. int i;
  807. double chisq;
  808. double last_chisq;
  809. double *a;
  810. double lambda;
  811. FILE *device;
  812. {
  813.     int k;
  814.  
  815.     fprintf(device, "\n\n\
  816. Iteration %d\n\
  817. WSSR        : %-15g   delta(WSSR)/WSSR   : %g\n\
  818. delta(WSSR) : %-15g   limit for stopping : %g\n\
  819. lambda      : %g\n\n%s parameter values\n\n",
  820.       i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq) / chisq : 0.0,
  821.         chisq - last_chisq, epsilon, lambda,
  822.         (i > 0 ? "resultant" : "initial set of free"));
  823.     for (k = 0; k < num_params; k++)
  824.     fprintf(device, "%-15.15s = %g\n", par_name[k], a[k]);
  825. }
  826.  
  827.  
  828.  
  829. /*****************************************************************
  830.     is_empty: check for valid string entries
  831. *****************************************************************/
  832. static TBOOLEAN is_empty(s)
  833. char *s;
  834. {
  835.     while (*s == ' ' || *s == '\t' || *s == '\n')
  836.     s++;
  837.     return (TBOOLEAN) (*s == '#' || *s == '\0');
  838. }
  839.  
  840.  
  841. /*****************************************************************
  842.     get next word of a multi-word string, advance pointer
  843. *****************************************************************/
  844. char *get_next_word(s, subst)
  845. char **s;
  846. char *subst;
  847. {
  848.     char *tmp = *s;
  849.  
  850.     while (*tmp == ' ' || *tmp == '\t' || *tmp == '=')
  851.     tmp++;
  852.     if (*tmp == '\n' || *tmp == '\0')    /* not found */
  853.     return NULL;
  854.     if ((*s = strpbrk(tmp, " =\t\n")) == NULL)
  855.     *s = tmp + strlen(tmp);
  856.     *subst = **s;
  857.     *(*s)++ = '\0';
  858.     return tmp;
  859. }
  860.  
  861.  
  862. /*****************************************************************
  863.     check for variable identifiers
  864. *****************************************************************/
  865. static TBOOLEAN is_variable(s)
  866. char *s;
  867. {
  868.     while (*s != '\0') {
  869.     if (!isalnum((int) *s) && *s != '_')
  870.         return FALSE;
  871.     s++;
  872.     }
  873.     return TRUE;
  874. }
  875.  
  876. /*****************************************************************
  877.     first time settings
  878. *****************************************************************/
  879. void init_fit()
  880. {
  881.     func.at = (struct at_type *) NULL;    /* need to parse 1 time */
  882. }
  883.  
  884.  
  885. /*****************************************************************
  886.         Set a GNUPLOT user-defined variable
  887. ******************************************************************/
  888.  
  889. void setvar(varname, data)
  890. char *varname;
  891. struct value data;
  892. {
  893.     register struct udvt_entry *udv_ptr = first_udv, *last = first_udv;
  894.  
  895.     /* check if it's already in the table... */
  896.  
  897.     while (udv_ptr) {
  898.     last = udv_ptr;
  899.     if (!strcmp(varname, udv_ptr->udv_name))
  900.         break;
  901.     udv_ptr = udv_ptr->next_udv;
  902.     }
  903.  
  904.     if (!udv_ptr) {        /* generate new entry */
  905.     last->next_udv = udv_ptr = (struct udvt_entry *)
  906.         gp_alloc((unsigned int) sizeof(struct udvt_entry), "fit setvar");
  907.     udv_ptr->next_udv = NULL;
  908.     }
  909.     safe_strncpy(udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name));
  910.     udv_ptr->udv_value = data;
  911.     udv_ptr->udv_undef = FALSE;
  912. }
  913.  
  914.  
  915.  
  916. /*****************************************************************
  917.     Read INTGR Variable value, return 0 if undefined or wrong type
  918. *****************************************************************/
  919. int getivar(varname)
  920. char *varname;
  921. {
  922.     register struct udvt_entry *udv_ptr = first_udv;
  923.  
  924.     while (udv_ptr) {
  925.     if (!strcmp(varname, udv_ptr->udv_name))
  926.         return udv_ptr->udv_value.type == INTGR
  927.         ? udv_ptr->udv_value.v.int_val    /* valid */
  928.         : 0;        /* wrong type */
  929.     udv_ptr = udv_ptr->next_udv;
  930.     }
  931.     return 0;            /* not in table */
  932. }
  933.  
  934.  
  935. /*****************************************************************
  936.     Read DOUBLE Variable value, return 0 if undefined or wrong type
  937.    I dont think it's a problem that it's an integer - div
  938. *****************************************************************/
  939. static double getdvar(varname)
  940. char *varname;
  941. {
  942.     register struct udvt_entry *udv_ptr = first_udv;
  943.  
  944.     for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
  945.     if (strcmp(varname, udv_ptr->udv_name) == 0)
  946.         return real(&(udv_ptr->udv_value));
  947.  
  948.     /* get here => not found */
  949.     return 0;
  950. }
  951.  
  952. /*****************************************************************
  953.    like getdvar, but
  954.    - convert it from integer to real if necessary
  955.    - create it with value INITIAL_VALUE if not found or undefined
  956. *****************************************************************/
  957. static double createdvar(varname, value)
  958. char *varname;
  959. double value;
  960. {
  961.     register struct udvt_entry *udv_ptr = first_udv;
  962.  
  963.     for (; udv_ptr; udv_ptr = udv_ptr->next_udv)
  964.     if (strcmp(varname, udv_ptr->udv_name) == 0) {
  965.         if (udv_ptr->udv_undef) {
  966.         udv_ptr->udv_undef = 0;
  967.         (void) Gcomplex(&udv_ptr->udv_value, value, 0.0);
  968.         } else if (udv_ptr->udv_value.type == INTGR) {
  969.         (void) Gcomplex(&udv_ptr->udv_value, (double) udv_ptr->udv_value.v.int_val, 0.0);
  970.         }
  971.         return real(&(udv_ptr->udv_value));
  972.     }
  973.     /* get here => not found */
  974.  
  975.     {
  976.     struct value tempval;
  977.     (void) Gcomplex(&tempval, value, 0.0);
  978.     setvar(varname, tempval);
  979.     }
  980.  
  981.     return value;
  982. }
  983.  
  984.  
  985. /*****************************************************************
  986.     Split Identifier into path and filename
  987. *****************************************************************/
  988. static void splitpath(s, p, f)
  989. char *s;
  990. char *p;
  991. char *f;
  992. {
  993.     register char *tmp;
  994.     tmp = s + strlen(s) - 1;
  995.     while (*tmp != '\\' && *tmp != '/' && *tmp != ':' && tmp - s >= 0)
  996.     tmp--;
  997.     strcpy(f, tmp + 1);
  998.     safe_strncpy(p, s, (size_t) (tmp - s + 1));
  999.     p[tmp - s + 1] = NUL;
  1000. }
  1001.  
  1002.  
  1003. /*****************************************************************
  1004.     write the actual parameters to start parameter file
  1005. *****************************************************************/
  1006. void update(pfile, npfile)
  1007. char *pfile, *npfile;
  1008. {
  1009.     char fnam[256], path[256], sstr[256], pname[64], tail[127], *s = sstr,
  1010.     *tmp, c;
  1011.     char ifilename[256], *ofilename;
  1012.     FILE *of, *nf;
  1013.     double pval;
  1014.  
  1015.  
  1016.     /* update pfile npfile:
  1017.        if npfile is a valid file name,
  1018.        take pfile as input file and
  1019.        npfile as output file
  1020.      */
  1021.     if (VALID_FILENAME(npfile)) {
  1022.     safe_strncpy(ifilename, pfile, sizeof(ifilename));
  1023.     ofilename = npfile;
  1024.     } else {
  1025. #ifdef BACKUP_FILESYSTEM
  1026.     /* filesystem will keep original as previous version */
  1027.     safe_strncpy(ifilename, pfile, sizeof(ifilename));
  1028. #else
  1029.     backup_file(ifilename, pfile);    /* will Eex if it fails */
  1030. #endif
  1031.     ofilename = pfile;
  1032.     }
  1033.  
  1034.     /* split into path and filename */
  1035.     splitpath(ifilename, path, fnam);
  1036.     if (!(of = fopen(ifilename, "r")))
  1037.     Eex2("parameter file %s could not be read", ifilename);
  1038.  
  1039.     if (!(nf = fopen(ofilename, "w")))
  1040.     Eex2("new parameter file %s could not be created", ofilename);
  1041.  
  1042.     while (fgets(s = sstr, sizeof(sstr), of) != NULL) {
  1043.  
  1044.     if (is_empty(s)) {
  1045.         fputs(s, nf);    /* preserve comments */
  1046.         continue;
  1047.     }
  1048.     if ((tmp = strchr(s, '#')) != NULL) {
  1049.         safe_strncpy(tail, tmp, sizeof(tail));
  1050.         *tmp = NUL;
  1051.     } else
  1052.         strcpy(tail, "\n");
  1053.  
  1054.     tmp = get_next_word(&s, &c);
  1055.     if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
  1056.         (void) fclose(nf);
  1057.         (void) fclose(of);
  1058.         Eex2("syntax error in parameter file %s", fnam);
  1059.     }
  1060.     safe_strncpy(pname, tmp, sizeof(pname));
  1061.     /* next must be '=' */
  1062.     if (c != '=') {
  1063.         tmp = strchr(s, '=');
  1064.         if (tmp == NULL) {
  1065.         (void) fclose(nf);
  1066.         (void) fclose(of);
  1067.         Eex2("syntax error in parameter file %s", fnam);
  1068.         }
  1069.         s = tmp + 1;
  1070.     }
  1071.     tmp = get_next_word(&s, &c);
  1072.     if (!sscanf(tmp, "%lf", &pval)) {
  1073.         (void) fclose(nf);
  1074.         (void) fclose(of);
  1075.         Eex2("syntax error in parameter file %s", fnam);
  1076.     }
  1077.     if ((tmp = get_next_word(&s, &c)) != NULL) {
  1078.         (void) fclose(nf);
  1079.         (void) fclose(of);
  1080.         Eex2("syntax error in parameter file %s", fnam);
  1081.     }
  1082.     /* now modify */
  1083.  
  1084.     if ((pval = getdvar(pname)) == 0)
  1085.         pval = (double) getivar(pname);
  1086.  
  1087.     sprintf(sstr, "%g", pval);
  1088.     if (!strchr(sstr, '.') && !strchr(sstr, 'e'))
  1089.         strcat(sstr, ".0");    /* assure CMPLX-type */
  1090.  
  1091.     fprintf(nf, "%-15.15s = %-15.15s   %s", pname, sstr, tail);
  1092.     }
  1093.  
  1094.     if (fclose(nf) || fclose(of))
  1095.     Eex("I/O error during update");
  1096.  
  1097. }
  1098.  
  1099.  
  1100. /*****************************************************************
  1101.     Backup a file by renaming it to something useful. Return
  1102.     the new name in tofile
  1103. *****************************************************************/
  1104. static void backup_file(tofile, fromfile)
  1105. char *tofile;
  1106. const char *fromfile;
  1107. /* tofile must point to a char array[] or allocated data. See update() */
  1108. {
  1109. #if defined (WIN32) || defined(MSDOS) || defined(VMS)
  1110.     char *tmpn;
  1111. #endif
  1112.  
  1113.     /* win32 needs to have two attempts at the rename, since it may
  1114.      * be running on win32s with msdos 8.3 names
  1115.      */
  1116.  
  1117. /* first attempt, for all o/s other than MSDOS */
  1118.  
  1119. #ifndef MSDOS
  1120.  
  1121.     strcpy(tofile, fromfile);
  1122.  
  1123. #ifdef VMS
  1124.     /* replace all dots with _, since we will be adding the only
  1125.      * dot allowed in VMS names
  1126.      */
  1127.     while ((tmpn = strchr(tofile, '.')) != NULL)
  1128.     *tmpn = '_';
  1129. #endif /*VMS */
  1130.  
  1131.     strcat(tofile, BACKUP_SUFFIX);
  1132.  
  1133.     if (rename(fromfile, tofile) == 0)
  1134.     return;            /* hurrah */
  1135. #endif
  1136.  
  1137.  
  1138. #if defined (WIN32) || defined(MSDOS)
  1139.  
  1140.     /* first attempt for msdos. Second attempt for win32s */
  1141.  
  1142.     /* beware : strncpy is very dangerous since it does not guarantee
  1143.      * to terminate the string. Copy up to 8 characters. If exactly
  1144.      * chars were copied, the string is not terminated. If the
  1145.      * source string was shorter than 8 chars, no harm is done
  1146.      * (here) by writing to offset 8.
  1147.      */
  1148.     safe_strncpy(tofile, fromfile, 8);
  1149.     /* tofile[8] = NUL; */
  1150.  
  1151.     while ((tmpn = strchr(tofile, '.')) != NULL)
  1152.     *tmpn = '_';
  1153.  
  1154.     strcat(tofile, BACKUP_SUFFIX);
  1155.  
  1156.     if (rename(fromfile, tofile) == 0)
  1157.     return;            /* success */
  1158.  
  1159. #endif /* win32 || msdos */
  1160.  
  1161.     /* get here => rename failed. */
  1162.     Eex3("Could not rename file %s to %s", fromfile, tofile);
  1163. }
  1164.  
  1165.  
  1166. /*****************************************************************
  1167.     Interface to the classic gnuplot-software
  1168. *****************************************************************/
  1169.  
  1170. void do_fit()
  1171. {
  1172.     TBOOLEAN autorange_x = 3, autorange_y = 3;    /* yes */
  1173.     /* HBB 980401: new: z range specification */
  1174.     TBOOLEAN autorange_z = 3;
  1175.     double min_x, max_x;    /* range to fit */
  1176.     double min_y, max_y;    /* range to fit */
  1177.     /* HBB 980401: new: z range specification */
  1178.     double min_z, max_z;    /* range to fit */
  1179.     int dummy_x = -1, dummy_y = -1;    /* eg  fit [u=...] [v=...] */
  1180.     /* HBB 981210: memorize position of possible third [ : ] spec: */
  1181.     int zrange_token = -1;
  1182.  
  1183.     int i;
  1184.     double v[4];
  1185.     double tmpd;
  1186.     time_t timer;
  1187.     int token1, token2, token3;
  1188.     char *tmp;
  1189.  
  1190.     /* first look for a restricted x fit range... */
  1191.  
  1192.     if (equals(c_token, "[")) {
  1193.     c_token++;
  1194.     if (isletter(c_token)) {
  1195.         if (equals(c_token + 1, "=")) {
  1196.         dummy_x = c_token;
  1197.         c_token += 2;
  1198.         }
  1199.         /* else parse it as a xmin expression */
  1200.     }
  1201.     autorange_x = load_range(FIRST_X_AXIS, &min_x, &max_x, autorange_x);
  1202.     if (!equals(c_token, "]"))
  1203.         int_error("']' expected", c_token);
  1204.     c_token++;
  1205.     }
  1206.     /* ... and y */
  1207.  
  1208.     if (equals(c_token, "[")) {
  1209.     c_token++;
  1210.     if (isletter(c_token)) {
  1211.         if (equals(c_token + 1, "=")) {
  1212.         dummy_y = c_token;
  1213.         c_token += 2;
  1214.         }
  1215.         /* else parse it as a ymin expression */
  1216.     }
  1217.     autorange_y = load_range(FIRST_Y_AXIS, &min_y, &max_y, autorange_y);
  1218.     if (!equals(c_token, "]"))
  1219.         int_error("']' expected", c_token);
  1220.     c_token++;
  1221.     }
  1222.     /* HBB 980401: new: allow restricting the z range as well */
  1223.     /* ... and z */
  1224.  
  1225.     if (equals(c_token, "[")) {
  1226.     zrange_token = c_token++;
  1227.     if (isletter(c_token))
  1228.         /* do *not* allow the z range being given with a variable name */
  1229.         int_error("Can't re-name dependent variable", c_token);
  1230.     autorange_z = load_range(FIRST_Z_AXIS, &min_z, &max_z, autorange_z);
  1231.     if (!equals(c_token, "]"))
  1232.         int_error("']' expected", c_token);
  1233.     c_token++;
  1234. #if 0 /* HBB 981210: move this to a later point */
  1235.     } else {
  1236.     /* Just in case I muck up things below: make sure that the z
  1237.      * range is the same as the y range, if it didn't get specified
  1238.      */
  1239.     autorange_z = autorange_y;
  1240.     min_z = min_y;
  1241.     max_z = max_y;
  1242. #endif
  1243.     }
  1244.  
  1245.  
  1246.     /* now compile the function */
  1247.  
  1248.     token1 = c_token;
  1249.  
  1250.     if (func.at) {
  1251.     free(func.at);
  1252.     func.at = NULL;        /* in case perm_at() does int_error */
  1253.     }
  1254.     dummy_func = &func;
  1255.  
  1256.  
  1257.     /* set the dummy variable names */
  1258.  
  1259.     if (dummy_x >= 0)
  1260.     copy_str(c_dummy_var[0], dummy_x, MAX_ID_LEN);
  1261.     else
  1262.     strcpy(c_dummy_var[0], dummy_var[0]);
  1263.  
  1264.     if (dummy_y >= 0)
  1265.     copy_str(c_dummy_var[0], dummy_y, MAX_ID_LEN);
  1266.     else
  1267.     strcpy(c_dummy_var[1], dummy_var[1]);
  1268.  
  1269.     func.at = perm_at();
  1270.     dummy_func = NULL;
  1271.  
  1272.     token2 = c_token;
  1273.  
  1274.     /* use datafile module to parse the datafile and qualifiers */
  1275.  
  1276.     columns = df_open(4);    /* up to 4 using specs allowed */
  1277.     if (columns == 1)
  1278.     int_error("Need 2 to 4 using specs", c_token);
  1279.  
  1280.     /* HBB 980401: if this is a single-variable fit, we shouldn't have
  1281.      * allowed a variable name specifier for 'y': */
  1282.     if ((dummy_y >= 0) && (columns < 4))
  1283.     int_error("Can't re-name 'y' in a one-variable fit", dummy_y);
  1284.  
  1285.     /* HBB 981210: two range specs mean different things, depending
  1286.      * on wether this is a 2D or 3D fit */
  1287.     if (columns<4) {
  1288.       if (zrange_token != -1)
  1289.     int_error("Three range-specs not allowed in on-variable fit", zrange_token);
  1290.       else {
  1291.     /* 2D fit, 2 ranges: second range is for *z*, not y: */
  1292.     autorange_z = autorange_y;
  1293.     min_z = min_y;
  1294.     max_z = max_y;
  1295.       }
  1296.     }
  1297.  
  1298.     /* defer actually reading the data until we have parsed the rest
  1299.      * of the line */
  1300.  
  1301.     token3 = c_token;
  1302.  
  1303.     tmpd = getdvar(FITLIMIT);    /* get epsilon if given explicitly */
  1304.     if (tmpd < 1.0 && tmpd > 0.0)
  1305.     epsilon = tmpd;
  1306.  
  1307.     /* HBB 970304: maxiter patch */
  1308.     maxiter = getivar(FITMAXITER);
  1309.  
  1310.     /* get startup value for lambda, if given */
  1311.     tmpd = getdvar(FITSTARTLAMBDA);
  1312.  
  1313.     if (tmpd > 0.0) {
  1314.     startup_lambda = tmpd;
  1315.     printf("Lambda Start value set: %g\n", startup_lambda);
  1316.     }
  1317.  
  1318.     /* get lambda up/down factor, if given */
  1319.     tmpd = getdvar(FITLAMBDAFACTOR);
  1320.  
  1321.     if (tmpd > 0.0) {
  1322.     lambda_up_factor = lambda_down_factor = tmpd;
  1323.     printf("Lambda scaling factors reset:  %g\n", lambda_up_factor);
  1324.     }
  1325.     *fit_script = NUL;
  1326.     if ((tmp = getenv(FITSCRIPT)) != NULL)
  1327.     strcpy(fit_script, tmp);
  1328.  
  1329.     tmp = getenv(GNUFITLOG);    /* open logfile */
  1330.     if (tmp != NULL) {
  1331.     char *tmp2 = &tmp[strlen(tmp) - 1];
  1332.     if (*tmp2 == '/' || *tmp2 == '\\')
  1333.         sprintf(logfile, "%s%s", tmp, logfile);
  1334.     else
  1335.         strcpy(logfile, tmp);
  1336.     }
  1337.     if (!log_f && /* div */ !(log_f = fopen(logfile, "a")))
  1338.     Eex2("could not open log-file %s", logfile);
  1339.     fputs("\n\n*******************************************************************************\n", log_f);
  1340.     (void) time(&timer);
  1341.     fprintf(log_f, "%s\n\n", ctime(&timer));
  1342.     {
  1343.     char *line = NULL;
  1344.  
  1345.     m_capture(&line, token2, token3 - 1);
  1346.     fprintf(log_f, "FIT:    data read from %s\n", line);
  1347.     free(line);
  1348.     }
  1349.  
  1350.     max_data = MAX_DATA;
  1351.     fit_x = vec(max_data);    /* start with max. value */
  1352.     fit_y = vec(max_data);
  1353.     fit_z = vec(max_data);
  1354.  
  1355.     /* first read in experimental data */
  1356.  
  1357.     err_data = vec(max_data);
  1358.     num_data = 0;
  1359.  
  1360.     while ((i = df_readline(v, 4)) != EOF) {
  1361.     if (num_data >= max_data) {
  1362.         /* increase max_data by factor of 1.5 */
  1363.         max_data = (max_data * 3) / 2;
  1364.         if (0
  1365.         || !redim_vec(&fit_x, max_data)
  1366.         || !redim_vec(&fit_y, max_data)
  1367.         || !redim_vec(&fit_z, max_data)
  1368.         || !redim_vec(&err_data, max_data)
  1369.         ) {
  1370.         /* Some of the reallocations went bad: */
  1371.         df_close();
  1372.         Eex2("Out of memory in fit: too many datapoints (%d)?", max_data);
  1373.         } else {
  1374.         /* Just so we know that the routine is at work: */
  1375.         fprintf(STANDARD, "Max. number of data points scaled up to: %d\n",
  1376.             max_data);
  1377.         }
  1378.     }
  1379.     switch (i) {
  1380.     case DF_UNDEFINED:
  1381.     case DF_FIRST_BLANK:
  1382.     case DF_SECOND_BLANK:
  1383.         continue;
  1384.     case 0:
  1385.         Eex2("bad data on line %d of datafile", df_line_number);
  1386.         break;
  1387.     case 1:        /* only z provided */
  1388.         v[2] = v[0];
  1389.         v[0] = (double) df_datum;
  1390.         break;
  1391.     case 2:        /* x and z */
  1392.         v[2] = v[1];
  1393.         break;
  1394.  
  1395.         /* only if the explicitly asked for 4 columns do we
  1396.          * do a 3d fit. (We can get here if they didn't
  1397.          * specify a using spec, and the file has 4 columns
  1398.          */
  1399.     case 4:        /* x, y, z, error */
  1400.         if (columns == 4)
  1401.         break;
  1402.         /* else fall through */
  1403.     case 3:        /* x, z, error */
  1404.         v[3] = v[2];    /* error */
  1405.         v[2] = v[1];    /* z */
  1406.         break;
  1407.  
  1408.     }
  1409.  
  1410.     /* skip this point if it is out of range */
  1411.     if (!(autorange_x & 1) && (v[0] < min_x))
  1412.         continue;
  1413.     if (!(autorange_x & 2) && (v[0] > max_x))
  1414.         continue;
  1415.     /* HBB 980401: yrange is only relevant for 3d-fits! */
  1416.     if (columns == 4) {
  1417.         if (!(autorange_y & 1) && (v[1] < min_y))
  1418.         continue;
  1419.         if (!(autorange_y & 2) && (v[1] > max_y))
  1420.         continue;
  1421.     }
  1422.     /* HBB 980401: check *z* range for all fits */
  1423.     if (!(autorange_z & 1) && (v[2] < min_z))
  1424.         continue;
  1425.     if (!(autorange_z & 2) && (v[2] > max_z))
  1426.         continue;
  1427.  
  1428.     fit_x[num_data] = v[0];
  1429.     fit_y[num_data] = v[1];
  1430.     fit_z[num_data] = v[2];
  1431.  
  1432.     /* we only use error if _explicitly_ asked for by a
  1433.      * using spec
  1434.      */
  1435.     err_data[num_data++] = (columns > 2) ? v[3] : 1;
  1436.     }
  1437.     df_close();
  1438.  
  1439.     if (num_data <= 1)
  1440.     Eex("No data to fit ");
  1441.  
  1442.     /* now resize fields to actual length: */
  1443.     redim_vec(&fit_x, num_data);
  1444.     redim_vec(&fit_y, num_data);
  1445.     redim_vec(&fit_z, num_data);
  1446.     redim_vec(&err_data, num_data);
  1447.  
  1448.     fprintf(log_f, "        #datapoints = %d\n", num_data);
  1449.  
  1450.     if (columns < 3)
  1451.     fputs("        residuals are weighted equally (unit weight)\n\n", log_f);
  1452.  
  1453.     {
  1454.     char *line = NULL;
  1455.  
  1456.     m_capture(&line, token1, token2 - 1);
  1457.     fprintf(log_f, "function used for fitting: %s\n", line);
  1458.     free(line);
  1459.     }
  1460.  
  1461.     /* read in parameters */
  1462.  
  1463.     max_params = MAX_PARAMS;    /* HBB 971023: make this resizeable */
  1464.  
  1465.     if (!equals(c_token, "via"))
  1466.     int_error("Need via and either parameter list or file", c_token);
  1467.  
  1468.     a = vec(max_params);
  1469.     par_name = (fixstr *) gp_alloc((max_params + 1) * sizeof(fixstr), "fit param");
  1470.     num_params = 0;
  1471.  
  1472.     if (isstring(++c_token)) {    /* It's a parameter *file* */
  1473.     TBOOLEAN fixed;
  1474.     double tmp_par;
  1475.     char c, *s;
  1476.     char sstr[MAX_LINE_LEN];
  1477.     FILE *f;
  1478.  
  1479.     quote_str(sstr, c_token++, MAX_LINE_LEN);
  1480.  
  1481.     /* get parameters and values out of file and ignore fixed ones */
  1482.  
  1483.     fprintf(log_f, "fitted parameters and initial values from file: %s\n\n", sstr);
  1484.     if (!(f = fopen(sstr, "r")))
  1485.         Eex2("could not read parameter-file %s", sstr);
  1486.     while (TRUE) {
  1487.         if (!fgets(s = sstr, (int) sizeof(sstr), f))    /* EOF found */
  1488.         break;
  1489.         if ((tmp = strstr(s, FIXED)) != NULL) {    /* ignore fixed params */
  1490.         *tmp = NUL;
  1491.         fprintf(log_f, "FIXED:  %s\n", s);
  1492.         fixed = TRUE;
  1493.         } else
  1494.         fixed = FALSE;
  1495.         if ((tmp = strchr(s, '#')) != NULL)
  1496.         *tmp = NUL;
  1497.         if (is_empty(s))
  1498.         continue;
  1499.         tmp = get_next_word(&s, &c);
  1500.         if (!is_variable(tmp) || strlen(tmp) > MAX_ID_LEN) {
  1501.         (void) fclose(f);
  1502.         Eex("syntax error in parameter file");
  1503.         }
  1504.         safe_strncpy(par_name[num_params], tmp, sizeof(fixstr));
  1505.         /* next must be '=' */
  1506.         if (c != '=') {
  1507.         tmp = strchr(s, '=');
  1508.         if (tmp == NULL) {
  1509.             (void) fclose(f);
  1510.             Eex("syntax error in parameter file");
  1511.         }
  1512.         s = tmp + 1;
  1513.         }
  1514.         tmp = get_next_word(&s, &c);
  1515.         if (sscanf(tmp, "%lf", &tmp_par) != 1) {
  1516.         (void) fclose(f);
  1517.         Eex("syntax error in parameter file");
  1518.         }
  1519.         /* make fixed params visible to GNUPLOT */
  1520.         if (fixed) {
  1521.         struct value tempval;
  1522.         (void) Gcomplex(&tempval, tmp_par, 0.0);
  1523.         /* use parname as temp */
  1524.         setvar(par_name[num_params], tempval);
  1525.         } else {
  1526.         if (num_params >= max_params) {
  1527.             max_params = (max_params * 3) / 2;
  1528.             if (0
  1529.             || !redim_vec(&a, max_params)
  1530.             || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
  1531.             ) {
  1532.             (void) fclose(f);
  1533.             Eex("Out of memory in fit: too many parameters?");
  1534.             }
  1535.         }
  1536.         a[num_params++] = tmp_par;
  1537.         }
  1538.  
  1539.         if ((tmp = get_next_word(&s, &c)) != NULL) {
  1540.         (void) fclose(f);
  1541.         Eex("syntax error in parameter file");
  1542.         }
  1543.     }
  1544.     (void) fclose(f);
  1545.  
  1546.     } else {
  1547.     /* not a string after via: it's a variable listing */
  1548.  
  1549.     fputs("fitted parameters initialized with current variable values\n\n", log_f);
  1550.     do {
  1551.         if (!isletter(c_token))
  1552.         Eex("no parameter specified");
  1553.         capture(par_name[num_params], c_token, c_token, (int) sizeof(par_name[0]));
  1554.         if (num_params >= max_params) {
  1555.         max_params = (max_params * 3) / 2;
  1556.         if (0
  1557.             || !redim_vec(&a, max_params)
  1558.             || !(par_name = gp_realloc(par_name, (max_params + 1) * sizeof(fixstr), "fit param resize"))
  1559.             ) {
  1560.             Eex("Out of memory in fit: too many parameters?");
  1561.         }
  1562.         }
  1563.         /* create variable if it doesn't exist */
  1564.         a[num_params] = createdvar(par_name[num_params], INITIAL_VALUE);
  1565.         ++num_params;
  1566.     } while (equals(++c_token, ",") && ++c_token);
  1567.     }
  1568.  
  1569.     redim_vec(&a, num_params);
  1570.     par_name = (fixstr *) gp_realloc(par_name, (num_params + 1) * sizeof(fixstr), "fit param");
  1571.  
  1572.     if (num_data < num_params)
  1573.     Eex("Number of data points smaller than number of parameters");
  1574.  
  1575.     /* avoid parameters being equal to zero */
  1576.     for (i = 0; i < num_params; i++)
  1577.     if (a[i] == 0)
  1578.         a[i] = NEARLY_ZERO;
  1579.  
  1580.     (void) regress(a);
  1581.  
  1582.     (void) fclose(log_f);
  1583.     log_f = NULL;
  1584.     free(fit_x);
  1585.     free(fit_y);
  1586.     free(fit_z);
  1587.     free(err_data);
  1588.     free(a);
  1589.     free(par_name);
  1590.     if (func.at) {
  1591.     free(func.at);        /* release perm. action table */
  1592.     func.at = (struct at_type *) NULL;
  1593.     }
  1594.     safe_strncpy(last_fit_command, input_line, sizeof(last_fit_command));
  1595. }
  1596.  
  1597. #if defined(ATARI) || defined(MTOS)
  1598. static int kbhit()
  1599. {
  1600.     fd_set rfds;
  1601.     struct timeval timeout;
  1602.  
  1603.     timeout.tv_sec = 0;
  1604.     timeout.tv_usec = 0;
  1605.     rfds = (fd_set) (1L << fileno(stdin));
  1606.     return ((select(0, &rfds, NULL, NULL, &timeout) > 0) ? 1 : 0);
  1607. }
  1608. #endif
  1609.